Comprehensive analysis of single cell and bulk data develops a promising prognostic signature for improving immunotherapy responses in ovarian cancer

The tumor heterogeneity is an important cause of clinical therapy failure and yields distinct prognosis in ovarian cancer (OV). Using the advantages of integrated single cell RNA sequencing (scRNA-seq) and bulk data to decode tumor heterogeneity remains largely unexplored. Four public datasets were enrolled in this study, including E-MTAB-8107, TCGA-OV, GSE63885, and GSE26193 cohorts. Random forest algorithm was employed to construct a multi-gene prognostic panel and further evaluated by receiver operator characteristic (ROC), calibration curve, and Cox regression. Subsequently, molecular characteristics were deciphered, and treatments strategies were explored to deliver precise therapy. The landscape of cell subpopulations and functional characteristics, as well as the dynamic of macrophage cells were detailly depicted at single cell level, and then screened prognostic candidate genes. Based on the expression of candidate genes, a stable and robust cell characterized gene associated prognosis signature (CCIS) was developed, which harbored excellent performance at prognosis assessment and patient stratification. The ROC and calibration curves, and Cox regression analysis elucidated CCIS could serve as serve as an independent factor for predicting prognosis. Moreover, a promising clinical tool nomogram was also constructed according to stage and CCIS. Through comprehensive investigations, patients in low-risk group were charactered by favorable prognosis, elevated genomic variations, higher immune cell infiltrations, and superior antigen presentation. For individualized treatment, patients in low-risk group were inclined to better immunotherapy responses. This study dissected tumor heterogeneity and afforded a promising prognostic signature, which was conducive to facilitating clinical outcomes for patients with OV.


Introduction
Ovarian cancer (OV) is a highly frequent malignant tumor among female populations [1].The mortality rate of OV is the highest among female gynecological cancers [2].Both high recurrence and easy drug tolerance are its classical characteristics, once diagnosed, current treatment choices often are limited [3,4].The systemic interventions of surgery, radiation therapy, and chemotherapy drugs (e.g., paclitaxel and carboplatin) are the mainstay of therapeutic strategies in clinical practice [5].However, due to tumor heterogeneity and unreasonable stratification, clinical treatments are still unsatisfactory, leading five-year survival rate less 50% [5].Notably, immunotherapy is gradually turn into emerging and burgeoning, which is also applied to numerous solid tumors, such as OV, hepatocellular carcinoma, lung adenocarcinoma [6,7].Whereas the application of immunotherapy is also failed to favorable clinical expectations in OV.Currently, effective stratification tool is still lacking, and patients under similar treatments perform distinct clinical outcomes, which hamper efficacy.Therefore, it is clinically crucial to develop biomarkers, identify "high-risk" patients, and provide treatment recommendations for OV patients.
In recent year, single cell RNA sequencing (scRNA-seq) as an emerging technology has obtain tremendous attention.The scRNA-seq cover the shortage in exploring the transcriptome heterogeneity as traditional RNA-seq only reflect the average expression in tissues [8].Based on this technology, researchers can perform transcriptomic analyses of individual cells, exploring tumor heterogeneity and cellular plasticity [8].Previous study had reported that inter-tumoral and intra-tumoral heterogeneity lead a major cause of clinical therapy failure [9].Meanwhile, high clinical and molecular heterogeneity of OV hamper personalized management, representing undertreatment or overtreatment [5,10].In clinical practice, treatment decisions mainly focus on pathological features and ignore molecular characteristics, such as genomic variation and tumor microenvironment (TME).Hereditary intra-tumoral heterogeneity impacts cell clonal evolution and genomic alterations drive distinct cancer evolutionary patterns, resulting different clinical outcomes [11].The TME plays fundamental roles in development and shape of tumor heterogeneity by regulating cell infiltration and secretory substance [12].To deepen understanding of heterogeneity, ever more machine learning approaches are used to tumor research.Moreover, these approaches have displayed conspicuous advantages at developing signatures.Currently, it is insufficient that using machine learning method integrally decode scRNA-seq and bulk data, classifying patients and further exploring tumor heterogeneity.Therefore, combining machine learning and scRNA-seq might a better choice to decipher tumor heterogeneity and seek optimal treatment strategies in OV.
In our study, the landscape of single cell subpopulation and functional characteristics were depicted.Using random forest (RSF) algorithm, an accuracy and robust cell characterized gene associated prognosis signature (CCIS) was constructed.Two external datasets were further enrolled to verify the excellent performance of CCIS by receiver operator characteristic (ROC), calibration curve, and Cox regression.Notably, the CCIS could serve as an independent factor for predicting overall survival (OS), progression-free survival (PFS), and diseasefree survival (DFS).To further decode tumor heterogeneity, our study also uncovered the molecular characteristics underlying CCIS, including potential biological function, genomics variations, antigen presentation, and immune infiltration.For personalized patient, the immunotherapy evaluation and drugs development were further performed, seeking optimal treatment strategies.Patients in low-risk group were suitable for immunotherapy and patients in high-risk group were sensitive to Pazopanib, Vinorelbine, Vorinostat, Shikonin, and FH535.Taken together, this work depicted tumor heterogeneity and offered a dependable prognostic signature for OV patients, which was conducive to enable clinical management and tailor treatment recommendation.

Data acquisition and processing
A single-cell RNA sequencing data and three independent transcriptome expression datasets were retrospectively enrolled in this study.The single-cell data, E-MTAB-8107 was accessed from ArrayExpress dataset (https://www.ebi.ac.uk/arrayexpress/).The corresponding expression profiles and clinical information were downloaded from TCGA portal and Gene Expression Omnibus (GEO) platform, including TCGA-OV, GSE63885, and GSE26193 cohorts.All patients were screened based on the inclusion and exclusion criteria: (a) primary OV; (b) survival time over 30 days; (c) No under any preoperative chemoradiotherapy.For genomics, somatic mutation data and copy number variation (CNV) information were obtained from TCGA website and FireBrowse dataset (http://firebrowse.org/),respectively.The CNV data was processed using GISTIC2.0 algorithm.To increase the homogeneity of expression data, the RNA-seq raw counts data were transformed to transcripts per kilobase million (TPM) format, and further converted to log-2 style.The GSE63885 and GSE26193 from GPL570 platform were processed and normalized via affy packages.
The immunotherapy sequencing data (GSE160755) was enrolled to assess the predictive performance of signature.In addition, three other immunotherapy cohorts were also retrieved into this study to further evaluate the immunotherapy efficacy in different groups, including GSE78220, GSE135222, and VanAllen cohorts.All these data were from GEO database and published article.The microarray expression profile was processed and normalized by affy and lumi packages according to different platforms.Based on the RECIST v1.1 standard, the responders and non-responders were defined by patients with complete response (CR)/partial response (PR) and patients with stable disease (SD) /progressive disease (PD), respectively.

Single cell annotation and type recognition
All single-cell data were originated from tumor tissue.The workflow in Seurat R toolkit (v4.0.6) was exploited to process single-cell data for dimension reduction, cells subpopulation clustering, and visualization [13,14].Among of these, some cells were excluded as the following criteria: (a) over 10% mitochondria genes; (b) fewer than 500 transcripts/cell.Then, the data was natural-log transformed and top 2000 variable genes were screened, which were served as input features for principal component analysis (PCA).Based on RunPCA function implemented in Seurat, PCA matrix harbored 50 components were used to dimensional reduction.Using the t-distributed Stochastic Neighbor Embedding (tSNE) algorithm, the data was further dimensional reduction, and clustering visualization was performed via RunTSNE function.The unique marker genes of various cell types were identified by FindAllMarkers function and classic cell markers were employed to annotate distinct cell types.For labeled cell types, specific expression genes were filtered by log2 FC> 1 and adjusted P-value< 0.05, and then marker genes were partial displayed on Vinplot, DotPlot, and Heatmap.The pseudo-time trajectory analysis of single cells was conducted by Monocle 2 package.The gene set variation analysis (GSVA) was applied to single-cell gene set enrichment analysis, exploring biological characteristics among distinct cell types.

The construction and validation of prognostic signature
The specific expression genes of different cell types were further deciphered at bulk transcriptome level.To better focus on clinical outcome, univariate Cox analysis was utilized to screen genes with potential prognostic significance, which possessed both P-value< 0.1 and all hazard ratio (HR)> 1 or < 1 across as least two datasets.The candidate genes were available in S1 Table .Based on the expression of these genes, a machine learning algorithm, random survival forest, was executed to develop a marker associated signature for OS.Two independently external cohorts were further used to validate the signature.In addition, the ROC and calibration curves were performed to elaborate the performance at predicting prognosis.

Prognosis evaluation and functional analysis
Apart from OS, the PFS and DFS were both common clinical endpoint events.The prognostic value of signature was further elucidated at PFS and DFS in TCGA-OV dataset.Multivariate Cox regression were used to assess whether risk score could serve as an independent prognostic factor for OV, when adjusted for usual clinical traits, such as age, stage, and grade.The heterogeneous clinical outcome was usually mirrored by distinct molecular characteristics and underlying biological behaviors.Patients were divided into high-risk and low-risk groups according to the median of score value.The fgsea package was employed to conduct gene set enrichment analysis (GSEA), identifying potential biological characteristics [15].Gene list was generated from genes ranked in descending order based on differential expression analysis, and gene sets were derived from go.v7.5.1 and kegg.v7.5.1 on MSigDB online.Using cluster-Profiler package, the top 10 pathways were presented, which harbored adjusted P-value< 0.05 and higher normalized enrichment score (NES).

The landscape of mutation and copy number variation
To depict the landscape of genomic alterations underlying two groups, we systemically explored somatic mutation and CNV data in TCGA-OV dataset.The mutation landscape with top 20 mutation frequencies were decoded and tumor mutation burden (TMB) of all samples was assessed, using maftools package.The amplification and deletion of CNV were further investigated by various indicators, including fraction of genome alteration (FGA), genomic gained (FGG), genome lost (FGL), arm gain/loss, and focal gain/loss.In addition, the GISTIC 2.0 pipeline also detailly deciphered altered regions based on gain/loss percentage and gistic score in high-risk and low-risk groups, respectively [16].

Immune infiltration and antigen presentation
Based on gene expression, the ssGSEA algorithm were employed to measure the abundances of 28 immune cells and further decipher immunological TME of each patient [17].Subsequently, the leukocyte fraction was exploited to reflect degree of inflammatory infiltration between high-risk and low-risk groups [18].Likewise, both homologous recombination defects (HRD) and neoantigens level were calculated and compared [18].The latent links with risk score were further dissected by correlation analysis.The expression of HLA molecules was compared, indicating the different capacity of antigen presentation.Additionally, the antigen presentation score (APS) and immunophenoscore (IPS) were further quantized to hint potential anti-tumor ability [17,19].

Immunotherapy assessment and drug development
To facilitate clinical benefits, we focus on immunotherapy and therapeutic drugs.The expression of co-stimulatory and co-inhibitory immune checkpoints (ICPs) was compared between two groups [20].T cell inflammatory signature (TIS) was a popular indicator for evaluating immunotherapy and the score was calculated via ssGSEA algorithm based on 18 inflammatory genes, which higher score performed better curative response [21].Additionally, two prevalent approaches were also implemented to estimate immunotherapeutic efficacy, encompassing TIDE and Submap [22,23].For identifying underlying therapeutic drugs, the pRRophetic package embedded linear ridge regression model, was utilized to measure half-maximal inhibitory concentration (IC50) [24].A lower IC50 value mapped on more sensitive to corresponding therapeutic drugs for OV patients.

Validation through quantitative real-time PCR (qRT-PCR)
The qRT-PCR was employed to quantify signature genes expression and further verify the accuracy of prognosis signature.Overall, 68 frozen tissues were obtained from OV patients in The First Affiliated Hospital of Zhengzhou University who underwent surgical therapy.All patients provided written informed consent and all these tissues were approved to use by the Ethics Committee Board of The First Affiliated Hospital of Zhengzhou University (2020-KY-365).Total RNA was isolated from OV tissues using RNAiso Plus reagent (Takara, Dalian, China) and RNA quality was estimated using a NanoDrop One C (Waltham, MA, USA).The complementary DNA (cDNA) was generated using 1μg of total RNA with reverse transcribed.Using SYBR Assay I Low ROX (Eurogentec, USA) and SYBR1 Green PCR Master Mix (Yeason, Shanghai, China), the qRT-PCR was conducted to evaluate gene expression of 13 signature genes.The ΔCT (Ct mRNA-Ct GAPDH) method was utilized to calculate expression value that was normalized to GAPDH, and further log2 transformed for follow-up analysis.The sequences of qRT-PCR primers were available in S2 Table.

Statistical analysis
The differential expression analysis was conducted by limma package.The Kaplan-Meier method and Log-rank test were exploited to evaluate the different OS, PFS, and DFS between high-risk and low-risk groups.Using timeROC and rms package, ROC and calibration curves were depicted to uncover the accuracy of prognostic signature.The correlations between two continuous variables were calculated and evaluated according to Spearman's correlation analysis.The survival package was applied in univariate and multivariate Cox regression analysis.The student's t-test/Wilcox rank-sum test was adopted to perform comparisons between two continuous variables.All data processing, statistical analysis, and plotting were accomplished in R software (v 4.1.2).All the source codes for each figure are available in S1 File.The adjusted P-value were corrected via Bonferroni method and P-value< 0.05 was considered statistically significant.

Cell clustering and typing of tumor tissues
The E-MTAB-8107 data was enrolled in this work, which offered tumor tissues and possessed high quality of single cell RNA sequencing information.The tSNE approach was performed to visualize the results of unsupervised dimensionality reduction and single cell subcluster analysis, and 27 cell clusters were classified (Fig 1A).Cell types were identified and recognized based on well-defined markers, such as C1QA, CAQB, and CD3D.These cells were classified into seven cell subpopulations, including Macrophage, Tumor cells, T cells, Fibroblasts Endothelial cells, B cells, and Dendritic cells (Fig 1B).To further favor the accuracy of cell annotation, the top two markers of these cell subpopulations were visualized as heatmap, bubble chart, and violin plot, which were consistent with previous studies (Fig 1C -1E).Meanwhile, Based on 50 Hallmark pathways, single cell enrichment analyses suggested these cell types broadly involved in cancer activity, especially Macrophage, Tumor cells, T cells, and Fibroblasts (Fig 2E).Single cell analysis substantiated these cell types might play crucial role in tumor progression.Thus, differentially expressed genes of main cell types were extracted by "FindAllMarkers" function, which were employed to subsequent analysis.

The development and validation of a robust prognostic signature
Using univariate Cox regression analysis, a total of 16 candidate genes were screened and then employed to model development (S1 Table ).Based on these gene expression, a robust signature for OS was developed by RSF machine learning approach.The relative importance value of model genes over 0.1 was reserved in our signature (S2A Fig), which was defined as CCIS.According to the median of score, patients were assigned into high-risk and low-risk groups.In the TCGA-OV cohort, Kaplan-Meier survival analysis suggested that patients in high-risk group performed dismal prognosis (P< 0.0001) (Fig 3A).Our CCIS were further validated in GSE63885 (P = 0.0270) and GSE26193 (P = 0.0001) cohorts, respectively (Fig 3B and 3C).The area under the ROC curves (AUC) and calibration curves were used to evaluate the accuracy and robustness of CCIS.Notably, the detailed AUCs with OS at 1/3/5 years were TCGA-OV (0.749/0.809/0.875),GSE63885 (0.846/0.704/0.710),and GSE26193 (0.666/0.779/0.731),respectively (Fig 3D -3F).The calibration curves with OS at 1/3/5 years were also demonstrated CCIS displayed excellent performance at predicting prognosis (Fig 3G -3I).
Likewise, patients in high-risk group also presented poor prognosis at both PFS and DFS (P< 0.0001) (Fig 4A and 4B).After adjusted for common clinical traits, multivariate Cox regression analysis substantiated the CCIS could serve as an independent factor for predicting OS, PFS, and DFS (Fig 4C -4E To further evaluate the robust of our signature, a clinical in-house cohort was employed and the risk score of each patient was calculated.The OV patients were classified into highrisk and low-risk groups according to the median of risk score.The Kaplan-Meier analysis of OS also suggested that the low-risk group harbored a favorable prognosis (S1A Fig) .Meanwhile, the AUCs for predicting OS at 1/3/5 years were in-house cohort (0.721/0.854/0.912)(S1B Fig) .The above results reconfirmed the excellent performance and proved the clinical value of our signature again.

The biological function and genomics alterations
Various biological characteristics existed links with heterogeneous clinical outcomes and crosstalk with TME [25].The underlying functional status of patients was further deciphered by GSEA analysis, which was conducive to decode tumor progression.Based on GO and KEGG pathways, patients in high-risk groups performed more malignant phenotype, such as ECM receptor interaction and pathways in cancer (Fig 5A).Meanwhile, the results illustrated immune and inflammatory associated process were mainly enriched in patients in low-risk groups, such as antigen processing and presentation (Fig 5A).In addition, genomics alterations with specific molecular characteristics also had fundamental influences on tumor progression [26].To depict the landscape of somatic mutation, top 20 high-frequency mutated genes were delineated by oncoplot (S2B Fig).Among of these, TP53, TTN, and CSMD3 were the first three and possessed tight relationships with tumorigenesis.However, there were no pronounced differences of mutation frequency between high-risk and low-risk groups  In addition, the CNV of patients were also dissected and evaluated between two groups.The genomic alterations were quantified in bases, fragments, and chromosome level.Interestingly, FGA, FGG, and FGL were compared and exhibited no difference between high-risk and low-risk groups (Fig 5F -5H).The genomics gain changes at chromosome arm were superior in high-risk group and genomics loss changes displayed similar tendency (Fig 6A and 6B).Both genomics gain and loss changes at chromosome focal were more conspicuous in low-risk group (Fig 6C -6D).Using GISTIC 2.0 pipeline, the gain/loss percentage and gistic score were all delineated, which highlighted the prominent differences of CNV at chromosome focal level (Figs 6E, 6F, S2C and S2D).Taken together, the CNV might bring heterogeneous biological patterns and clinical outcomes between two groups.

The performance of immune infiltration and antigen presentation
There were prominently distinct molecular characteristics and clinical outcomes in two groups, such as biological pathways, CNV, and prognosis.To yield knowledge about TME features, we explored the immune cell infiltration and evaluated tumor antigen presentation, further proving theoretical foundation for clinical treatments.As immune cells were key composition of TME and performed antitumor killing ability, the abundance of 28 immune cell infiltrations were firstly assessed by ssGSEA algorithm (Fig 7A).Major antitumor effector cells were elevated abundance in the low-risk group, including activated CD4 T cell and CD8 T cell (Fig 7A).In parallel, the leukocyte fraction was also superior in low-risk group relative high-risk group, hinting at more enriched inflammatory reserves (Fig 7B).Correlation analysis suggested negative relationships between risk score and HRD as well as neoantigens, respectively (Fig 7C and 7D).The higher HRD was eligible for more enormous neoantigens, which was likely to provide more backup resources for antitumor activity.
Based on HLA molecules, patients in low-risk group displayed elevated expression, such as HLA-A, HLA-B, HLA-C, and HLA-DRA (Fig 7E).Moreover, the APS and IPS were measured and compared in two groups.Consistent with the expression of HLA molecules, patients in low-risk group performed higher APS, indicating more powerful antigen presentation ability (Fig 7F).In addition, patients in low-risk group also presented superior IPS, suggesting higher immunogenicity and more sensitive to immunotherapy (Fig 7G ).

The exploration of treatments strategies to facilitate clinical management
In clinical practice, the ultimate goal was aiming to deliver curative treatments and improve prognosis [27].To bridge this gap, we further explored immunotherapy and therapeutic drugs to propose new insights on feasible treatments strategies.Apart from immune infiltration and antigen presentation, our study extended investigations on co-stimulatory and co-inhibitory molecules (Figs 8A and S3A).Our results substantiated both immune co-stimulatory and coinhibitory molecules were wealthier in low-risk group, such as CD70, ICOS, TNFSF13B, CD48 and VTCN1 (Figs 8A and S3A).Patients in low-risk group presented higher TIS score, implying better immunotherapy benefits (Fig 8B).Subsequently, the TIDE and Submap analysis were used to further favor populations with favorable immunotherapy response.The lower TIDE score means the better clinical efficacy, thus patients in low-risk groups were more compatible for immunotherapy (Fig 8C).Based on the similarity of expression matrix, SubMap analysis also indicated patients in low-risk group harbored more sensitive to anti-PD1 therapy (Fig 8D).Collectively, immunotherapy was recommended in low-risk group.The GZMA, GZMB, and CD107A represent the cytotoxic phenotype of immune cells, such as NK cells and T cells.We also used these molecules to evaluate the immunotherapy effect in OV patients.The expression of these molecules represented cytotoxic phenotype was higher in low-risk group (S4A-S4C Fig) .Thus, patients in low-risk group performed better responses on immunotherapy in GSE160755.For further evaluating the predictive performance of our signature, three other immunotherapy cohorts were enrolled, including GSE78220, GSE135222, and VanAllen cohorts.Notably, patients in the low-risk group exhibited a superior response to immunotherapy in all cohorts (S4D-S4F Fig) .The area under the curve values for predicting the accuracy of immunotherapeutic efficacy were 0.682/0.645/0.674 in GSE78220, GSE135222, and VanAllen cohorts, respectively (S4G-S4I Fig) .In conclusion, the results substantiated the risk scoring system was a promising immunotherapy indicator in OV patients.To tailor therapy strategies for high-risk group, the ridge regression model was exploited to develop underlying therapeutic drugs.Using pRRophetic package, the IC50 value of each patient was calculated and estimated, which exhibited lower value standing for more effective.Our results revealed that patients in high-risk group displayed superior sensitive to Pazopanib, Vinorelbine, Vorinostat, Shikonin, and FH535 (Fig 8E -8I).Patients in low-risk group were warranted that harnessing elevated response to Gefitinib, GNF.2, and Mitomycin.C (S3B-S3D Fig) .Overall, treatments strategies were provided for different populations, which might bring desirable efficacy for OV patients.

Discussion
Accumulating evidence suggests that inter-tumoral and intra-tumoral heterogeneity are a crucial cause of clinical therapy failure and yield distinct prognosis.The combinative explorations of scRNA-seq data and bulk data harbor tremendous advantages in decoding tumor heterogeneity, developing prognostic signature, making patients stratification, and providing treatment recommendation.To some extent Current used clinical management tool ignores molecular characteristics which hampers precise treatment, leading undertreatment or overtreatment [28].Therefore, developing a novel multi genes signature to explore tumor heterogeneity and seek optimal treatment strategies is imperative.In the present study, the landscape of cell subpopulations and functional characteristic was depicted.Then, the differentially expressed genes of main cell types were identified and further screened for signature discovery.Random forest is a popular machine learning approach and contains multiple decision trees and emphasizes important eigenvalue, which ensures signature with robust performance and strong extrapolation.Using RSF algorithm, a total of 13 genes were enrolled in our cell characterized gene associated prognosis signature, including WASF2, UBB, NR4A3, FBLN2, SH3PXD2A, MFAP4, UQCC2, QSOX1, LUZP1, SDF2L1, SFRP2, PRRX1, GLT8D2.Interestingly, all genes within CCIS have been demonstrated that tight correlations exist with tumor progression and prognosis across diverse cancer types.For instance, WASF2 is reported to be closely associated with TME and microsatellite instability, and elevated expression with poor prognosis in OV [29].The knockdown of UBB triggers cell apoptosis and β-catenin decreasing, which could serve as pro-survival gene in gastric adenocarcinoma [30].In addition, GSE63885 and GSE26193 cohorts further validated the robustness and accuracy of CCIS.Meanwhile, the CCIS could serve as an independent factor for predicting OS, PFS, and DFS.Therefore, our signature could work as a promising clinical tool for evaluating prognosis in OV patients.
Ever more studies have discovered that tumor heterogeneity could be mirrored by distinct molecular characteristics [31].Our study classified patients into high-risk group and low-risk group to further decode reciprocal correlations between tumor heterogeneity and molecular characteristics.As somatic mutation was similar, but CNV was mainly occurred in low-risk group, especially chromosome focal level, the CNV might be an important clue on for resulting distinct prognosis.Meanwhile, CNV was driving classical phenotypic variation and tumor progression via modifying and rearranging genes [32].Previous had also reported CNV affects the activity of various cancer-related pathways and favors tumor development [33].For example, the CNV amplification of PIK3CA enhanced PI3-kinase activity and elevated phosphorylated Akt level, leading aberrant cell proliferation and poor prognosis in gastric cancer [34].In addition, patients in low-risk group harbored higher HRD score, which might gain more clinical benefits from PARPi therapy [35].Elevated TMB was also the characteristic of low-risk group, harnessing more potential to immunotherapy response.To identify optimal beneficial populations, it is imperative to further explore the treatments strategies.
Our study further deciphered the performance of immune infiltration and antigen presentation.Patients in low-risk group was tended to harbor strong antitumor ability and possessed excellent resource for immune checkpoint inhibitor (ICI) treatment due to stored numerous immune killing cells, such as activated CD4 T cell and CD8 T cell.The IPS maps onto tumor immunogenicity while leukocyte fraction facilitates cytolytic activity, and both are favorable indicators for evaluating immunotherapy.Patients in low-risk group displayed higher IPS and leukocyte fraction, hinting better immunotherapeutic response.Previous study had reported that HRD triggers the generation of neoantigen epitopes and strengthens antitumor activity [36].Neoantigen-based therapies are gradually opening a new era in tumor immunotherapy [37].As expected, patients in low-risk group with superior tumor neoantigens indicates immunotherapy might be a worth considering approach for these populations.It is worth noting that the evaluation of immunotherapy efficacy is a complex process and should consider multiple steps [38].Patients in low-risk group presented higher expression of co-stimulatory and co-inhibitory ICPs, such as CD70, ICOS, TNFSF13B, CD48 and VTCN1.Among of them, CD70 as a novel immunotherapy target, promotes T cell and B cell activation and regulates immune response through interaction with CD27 [39].Apart from ICPs, antigen presentation is also a crucial step, which is broadly demonstrated to strengthen antitumor capacity [40].The superior APS and elevated expression of HLA genes also indicated patients in low-risk group adopted immunotherapy might arise in more curative benefits.Subsequently, three popular algorithms, TIS, TIDE, and Submap, were employed to further verify the above assumption.Consistently, all results declared patients in low-risk group could obtain thrilling curative benefits from immunotherapy.
Although our CCIS displays excellent performance at prognosis assessment and patient stratification, how to implement curative treatment for high-risk group with dismal prognosis is still an imperative issue.Compared to low-risk group, the immunotherapy efficacy was unfavorable in high-risk group owing to lower immune cells abundance and inferior antigen presentation.Thus, it might be a great option that develops potential drugs for high-risk group using ridge regression model.Based on expression data and drug sensitivity information, five latent sensitive drugs were determined for high-risk group, encompassing Pazopanib, Vinorelbine, Vorinostat, Shikonin, and FH535.For OV patients, the Pazopanib is one multitarget tyrosine inhibitor that could significantly prolong progression-free survival by combined utilization with paclitaxel [41].An experimental study suggests that the combination of Vinorelbine and ATM-3507 targeted Tpm3.1 compound could cause cell apoptosis and a hopeful treatment strategy for OV [42].Taken together, the above potential therapeutic drugs throw novel insights on treatment recommendations in high-risk group, facilitating clinical management.
Based on scRNA-seq and bulk data, a robust prognostic signature, CCIS was ultimately constructed to decode tumor heterogeneity.Although our CCIS harbored pronounced advantages at prognosis assessment and patient stratification, some limitations should be clarification.(a) all samples enrolled in this study were retrospective, thus future validation of CCIS should be conducted in multi-center prospective studies; (b) owing to lacked eligible patients with ICIs information, much more immunotherapy clinical research for OV patients should be retrieved to further favor our conclusions; (c) despite potential sensitive drugs were identified, their mode of actions remains further experiments.

Conclusions
This study deciphered tumor heterogeneity and established a feasible and robust prognostic signature in OV.With the delineation and exploration of molecular characteristics underlying the signature, CCIS, the treatment recommendations were provided for personalized population.Patients in high-risk group were more likely to gain clinical from benefits from potential therapeutic drugs, and patients in low-risk group were tailored immunotherapy that might

Fig 1 .
Fig 1. Cell clusters and cell subpopulation of ovarian cancer tissues.(A-B) The t-distributed stochastic neighbor embedding (tSNE) plot showing the cell clusters (A) and cell annotations (B) for distinct cell types in ovarian cancer.(C-E) The heatmap (C), bubble chart (D), and Violin plots (E) demonstrating the identity of each subpopulation through showing the expression of each cell type specific markers.https://doi.org/10.1371/journal.pone.0298125.g001 ).Overall, all results shed light on CCIS might be a promising clinical tool for patient stratification and prognosis assessment in OV patients.Considering the potential clinical transformation of CCIS, a stable nomogram was developed to serve clinical practice based on two indicators included stage and risk group (Fig 4F).The combination of clinical stage and CCIS exhibited superior accuracy and stability, which harbored AUCs with OS at 1/3/5 years were 0.781/0.820/0.875(Fig 4G).

Fig 2 .
Fig 2. The dynamic change of macrophage cells and potential biological process in single cell level.(A) Proportions of seven cell types originated from each ovarian cancer sample.(B) The t-SNE plot of only macrophage cells, displaying the cell clusters in ovarian cancer.(C-D) Differentiation trajectory of macrophage cells in HCC, with a color code for cell subpopulations (C) and pseudo-time (D).(E) Heatmap of 50 Hallmark gene sets among seven cell subpopulations using the gene set variation analysis (GSVA) algorithm.https://doi.org/10.1371/journal.pone.0298125.g002

Fig 4 .
Fig 4. The potential significance of clinical transformation.(A) Kaplan-Meier curves of progression-free survival (PFS) according to the signature in TCGA-OV cohort.(B) Kaplan-Meier curves of disease-free survival (DFS) according to the signature in TCGA-OV cohort.(C-E) Multivariate Cox regression of OS (C), PFS (D), and DFS (E) in TCGA-OV cohort.(F) A prognostic nomogram based on two indicators included stage and risk group.(G) Based on the nomogram, the Time-dependent ROC analysis for predicting OS at 1/3/5 years.https://doi.org/10.1371/journal.pone.0298125.g004

Fig 5 .
Fig 5.The underlying biological behaviors and genomic variations.(A) Based on GO and KEGG, exploring underlying biological characteristics by GSEA analysis according to the signature.(B-D) Composition mutated percentage of top three common somatic mutation gene, including TP53 (B), TTN (C), and CSMD3 (D) genes.(E) Distribution of tumor mutational burden (TMB) between high-risk and low-risk groups.(F-H) Distributions of fraction of genome alteration (F), fraction of genomic gained (G), and fraction of genome lost (H) between high-risk and low-risk groups.ns P> 0.05, **P< 0.01.https://doi.org/10.1371/journal.pone.0298125.g005

Fig 6 .
Fig 6.The landscape of copy number variation in ovarian cancer patients.(A-D) Distributions of copy number variations at chromosome arm gain (A), arm loss (B), focal gain (C), and focal loss (D) level, respectively.(E) The gain and loss of copy number percentage in high-risk group.(F) The gain and loss of copy number percentage in low-risk group.ns P> 0.05, *P< 0.05, ***P< 0.001.https://doi.org/10.1371/journal.pone.0298125.g006

Fig 7 .
Fig 7. Evaluation of immune infiltration and antigen presentation.(A) Infiltration abundance of 28 immune cell subsets (B) Comparison of leukocyte fraction between high-risk and low-risk groups.(C) The correlation between homologous recombination defects (HRD) and risk score, and the comparison of HRD between high-risk and low-risk groups.(D) The correlation between neoantigens and risk score, and the comparison of neoantigens between high-risk and low-risk groups.(E) Distribution of nine HLA molecular expressions.(F-G) Distribution of antigen presentation score (F) and immunophenoscore (G) between high-risk and low-risk groups.*P< 0.05, **P< 0.01, ***P< 0.001, ****P< 0.0001.https://doi.org/10.1371/journal.pone.0298125.g007